######## 
########  Load data again (you can start here if you don't want to reconstruct the basic files)
######## 

library(AER)
library(nleqslv)
library(readstata13)
library(matrixStats)
library(estimatr)

v2<-read.csv("Data/Temp/Farmer_Replication.csv",header=TRUE)	

y<-read.csv("Data/Temp/Agrovet_Replication.csv",header=TRUE)
z<-read.csv("Data/Temp/Costs_Shares_Replication.csv",header=TRUE)	
w<-read.csv("Data/Temp/Costs_Shares_Replication.csv",header=TRUE)	

a2<-as.matrix(read.csv("Data/Temp/Adoption_Replication.csv",header=TRUE))	
f2<-as.matrix(read.csv("Data/Temp/Expenditure_Replication.csv",header=TRUE))		
id<-as.matrix(read.csv("Data/Temp/ID_Replication.csv",header=TRUE))		

av_village<-as.character(tapply(as.character(v2$village_name),paste(v2$village,v2$ward,v2$district),FUN=unique))
av_ward<-as.character(tapply(as.character(v2$ward),paste(v2$village,v2$ward,v2$district),FUN=unique))
av_district<-as.character(tapply(as.character(v2$district),paste(v2$village,v2$ward,v2$district),FUN=unique))
av_market<-as.character(tapply(as.character(v2$market),paste(v2$village,v2$ward,v2$district),FUN=unique))
av_pop_market<-as.numeric(tapply(v2$census_bothsex,paste(v2$village,v2$ward,v2$district),FUN=sum,na.rm=TRUE))
popframe<-data.frame(av_village,av_ward,av_district,av_market,av_pop_market)

#Merge populations to agrovet data

y<-merge(y,popframe,all.x=T)

v2$census_bothsex<-v2$census_bothsex/sum(v2$census_bothsex,na.rm=TRUE)

z$same<-ifelse(z$survey_region_O==z$survey_region_D&z$survey_district_O==z$survey_district_D&z$survey_ward_O==z$survey_ward_D&z$survey_village_O==z$survey_village_D,1,0)


# Calculate rural and main-road distances

z$DIST_km_direct_rural<-z$DIST_km_direct*z$rural_km_share	
z$DIST_km_direct_main<-z$DIST_km_direct*(1-z$rural_km_share)

#Ensure factors are removed from variables

v2$region<-as.character(v2$region)
v2$district<-as.character(v2$district)
v2$market<-as.character(v2$market)
v2$ward<-as.character(v2$ward)
v2$village_name<-as.character(v2$village_name)

y$av_region<-as.character(y$av_region)
y$av_district<-as.character(y$av_district)
y$av_ward<-as.character(y$av_ward)
y$av_village<-as.character(y$av_village)

for(i in 1:10){
  z[,i]<-as.character(z[,i])
}

###  Construct bilateral matrix of distance.  J (rows) will be agrovets, and I (columns) will be villages

J<-nrow(y)
I<-nrow(v2)

d<-matrix(nrow=J,ncol=I) # distance matrix
r<-matrix(nrow=J,ncol=I) # rural share matrix


#construct distance matrices

for(i in 1:I){
  
  sz<-subset(z,survey_region_O==v2$region[i]&survey_district_O==v2$district[i]&survey_ward_O==v2$ward[i]&survey_village_O==v2$village[i])
  
  for(j in 1:J){	
    
    ssz<-subset(sz,survey_region_D==toupper(y$av_region[j])&survey_district_D==y$av_district[j]&survey_ward_D==y$av_ward[j]&survey_village_D==y$av_village[j])
    
    if(nrow(ssz)>0){
      d[j,i]<-ssz$DIST_km_direct[1]
      r[j,i]<-ssz$rural_km_share[1]				
    }
    if(nrow(ssz)==0){
      d[j,i]<-NA
      r[j,i]<-NA				
    }	
  }
  if(i%%100==0){print(i)}
}

minmat<-matrix(ncol=1,nrow=ncol(d))

for(i in 1:ncol(d)){
  minmat[i,1]<-min(d[,i],na.rm=TRUE)
}

rowrem<-grep(1,ifelse(rowSums(is.na(d))>500,1,0))  #agrovets with incomplete village linkages - 500 is just an arbitrarily large number chosen based on inspection of the data
colrem<-grep(1,ifelse(colSums(is.na(d))>300,1,0))  #villages with incomplete agrovet linkages - 300 is just an arbitrarily large number chosen based on inspection of the data

d2<-d[-rowrem,-colrem]  	#remove incomplete villages and agrovets
r2<-r[-rowrem,-colrem]

# farmer level adoption, expenditures, and IDs for counterfactuals

a2<-a2[,-colrem]
f2<-f2[,-colrem]
id<-id[,-colrem]

adopt_mat<-a2
omega_mat<-f2  #assigned initially as expenditures, but make into expenditure shares

for(i in 1:nrow(omega_mat)){
  omega_mat[i,]<-omega_mat[i,]/colSums(f2,na.rm=TRUE)
}

#Caclulate aggregate expenditures at the village level by taking the average within the sample and multiplying by village size

E_fert<-v2$expense_avg*v2$census_bothsex

#Agrovet fertilizer revenues

V_fert<-y$av_fert_rev_USD_w1

#define revenues as row vectors, and remove the agrovets for which there is no distance data for their village

V_fert<-as.matrix(V_fert[-rowrem])

#define adoption and expenditures as column vector.  Adoption of fertilizer includes with and without seeds.

E_fert<-t(as.matrix(E_fert[-colrem]))
E_fert<-ifelse(is.na(E_fert),0,E_fert)

####################            ######################## 
#### Calculate estimated trade costs for fertilizer ####
####################            ########################  

#Load data from MNL
xtau<-read.dta13("Data/buying_fert_MNL_binned_replication_edit.dta")
#xtau<-read.dta13("/Users/aspearot/Dropbox/MaLT 2017/Data/MaLT Project Data/Final Analysis Datasets/buying_fert_MNL_binned_032822.dta")

#calculate matrices of main road shares and rural road shares
m2<-d2*(1-r2)
s2<-d2*(r2)

#calculate elasticity adjusted trade costs estimates
mtau<-ifelse(m2,NA,NA)
mtau<-ifelse(m2>75,exp(as.numeric(xtau[9])),mtau)
mtau<-ifelse(m2>50&m2<=75,exp(as.numeric(xtau[8])),mtau)
mtau<-ifelse(m2>40&m2<=50,exp(as.numeric(xtau[7])),mtau)
mtau<-ifelse(m2>30&m2<=40,exp(as.numeric(xtau[6])),mtau)
mtau<-ifelse(m2>20&m2<=30,exp(as.numeric(xtau[5])),mtau)
mtau<-ifelse(m2>15&m2<=20,exp(as.numeric(xtau[4])),mtau)
mtau<-ifelse(m2>10&m2<=15,exp(as.numeric(xtau[3])),mtau)
mtau<-ifelse(m2>5&m2<=10,exp(as.numeric(xtau[2])),mtau)
mtau<-ifelse(m2>0&m2<=5,exp(as.numeric(xtau[1])),mtau)
mtau<-ifelse(m2<=0,1,mtau)

rtau<-ifelse(s2,NA,NA)
rtau<-ifelse(s2>20,exp(as.numeric(xtau[14])),rtau)
rtau<-ifelse(s2>15&s2<=20,exp(as.numeric(xtau[13])),rtau)
rtau<-ifelse(s2>10&s2<=15,exp(as.numeric(xtau[12])),rtau)
rtau<-ifelse(s2>5&s2<=10,exp(as.numeric(xtau[11])),rtau)
rtau<-ifelse(s2>0&s2<=5,exp(as.numeric(xtau[10])),rtau)
rtau<-ifelse(s2<=0,1,rtau)	

tau<-mtau*rtau

#### Calibration Step 1 - Normalize Total fertilizer expenditures and total fertilizer revenues to one.

E_fert2<-E_fert/sum(E_fert,na.rm=TRUE)
V_fert2<-V_fert/sum(V_fert,na.rm=TRUE)

#Define matrix for use in the contraction mapping

Ebig<-matrix(NA,nrow=nrow(d2),ncol=ncol(d2))

for(j in 1:nrow(Ebig)){
  Ebig[j,]<-E_fert2
}

Dpos<-as.matrix(ifelse(c(V_fert2)==0,0,1))

contract_start<-function(d){
  d2<-Dpos*as.matrix(d)
  dmatpre1<-d2
  
  dmat1<-matrix(dmatpre1,nrow=nrow(tau),ncol=ncol(tau))	
  #define lambda for agrovet j but without agrovet j in the numerator
  mshares1<-t(t((tau))/colSums(dmat1*(tau)))
  
  #predict the delta value for agrovet j
  dnew1<-(V_fert2/(rowSums(mshares1*Ebig)))
  
  #apply the normalization of deltas	
  dnew1<-dnew1/sum(dnew1)
  
  return(dnew1)
}

init_d0<-V_fert2
init_d<-V_fert2

for(i in 1:5000){
  
  new_d<-contract_start(init_d)
  if(i%%1==0){print(c(i,round(max(abs(new_d-init_d),na.rm=TRUE),digits=3)*1000000,cor(init_d,init_d0),max(init_d)))}
  if(max(abs(new_d-init_d),na.rm=TRUE)<1e-16){
    print(c(i,round(max(abs(new_d-init_d),na.rm=TRUE),digits=3)*1000000,cor(init_d,init_d0),max(init_d)))
    break}
  init_d<-new_d
}

share_func<-function(d){
  #d<-init_d1
  d2<-as.matrix(d)
  dmatpre1<-d2
  
  dmatpre1<-dmatpre1/sum(dmatpre1)
  
  dmat1<-matrix(dmatpre1,nrow=nrow(tau),ncol=ncol(tau))
  
  m1<-dmat1*(tau)
  mshares1<-t(t(m1)/colSums(m1))
  
  X1<-(V_fert2-rowSums(mshares1*Ebig))
  
  #SSR<-sum(X3^2,na.rm=TRUE)
  
  return(X1)		
}

#confirm solution by showing that the distance between revenues and expected expenditures at each agrovet is very small at the solution

max(abs(share_func(init_d)))

dmatpre1<-as.matrix(init_d)
dmatpre1<-dmatpre1/sum(dmatpre1)
dmat1<-matrix(dmatpre1,nrow=nrow(tau),ncol=ncol(tau))

m1<-dmat1*(tau)
mshares1<-t(t(m1)/colSums(m1))

Phi_f<-colSums(m1)

v3<-v2[-colrem,]
v3$farmers<-16

v3$Phi_f<-as.numeric(Phi_f)
v3$google_vil_city_km_std<-(v3$google_vil_city_km-mean(v3$google_vil_city_km,na.rm=TRUE))/sqrt(var(v3$google_vil_city_km,na.rm=TRUE))

h<-read.dta13("Data/Village Distances.dta")

h$market<-NULL

v3$ord<-seq(1,nrow(v3),1)
v4<-merge(v3,h,all.x=T)
v4<-v4[order(v4$ord),]

y2<-y[-rowrem,]

y2$deltas<-as.numeric(init_d)	#Quality adjusted prices, "delta", from the paper

#### RUN IV REGRESSIONS

elasreg0<-lm(log(deltas)~log(av_urea_rprice_USD_w1)+log(exper)+av_district,y2)		
elas_reduced<-(-as.numeric(coef(elasreg0)[2]))
elas_reduced

elasreg4<-ivreg(log(deltas)~log(av_urea_rprice_USD_w1)+log(exper)+av_district|log(av_urea_wprice_USD_w1)+log(exper)+av_district,data=y2)
elas_old<-(-as.numeric(coef(elasreg4)[2]))
elas_old

elasreg_IV<-ivreg(log(deltas)~log(av_urea_rprice_USD_w1)+log(exper)+av_district|log(av_urea_wprice_15_USD_w1)+log(exper)+av_district,data=subset(y2,av_pop_market>0))
elas_IV<-(-as.numeric(coef(elasreg_IV)[2]))
elas_IV

elasreg_IV_old<-ivreg(log(deltas)~log(av_urea_rprice_USD_w1)+log(exper)+av_district|log(av_pop_market)+log(exper)+av_district,data=subset(y2,av_pop_market>0))
elas_IV_old<-(-as.numeric(coef(elasreg_IV_old)[2]))
elas_IV_old		

#save for future use

write.csv(y2,"Results/IVDataset_Replication_Update.csv",row.names=FALSE,na=".")

zero<-function(l){
  ifelse(is.na(l),0,l)
}

infna<-function(l){
  ifelse(l==Inf,NA,l)
}

deltas_before<-y2$deltas
price_before<-y2$av_urea_rprice_USD_w1

#define elasticity used for analysis
elas<-elas_IV

#calculate main counterfactual change in trade costs
iceberg0<-tau^(-1/elas)-1
iceberg_half<-iceberg0/2
tau_half<-(1+iceberg_half)^(-elas)

#calculate counterfactual change in trade costs by main and rural roads separately
miceberg0<-mtau^(-1/elas)-1
miceberg_half<-miceberg0/2
riceberg0<-rtau^(-1/elas)-1
riceberg_half<-riceberg0/2	
tau_mainhalf<-(1+miceberg_half)^(-elas)*(1+riceberg0)^(-elas)	
tau_ruralhalf<-(1+miceberg0)^(-elas)*(1+riceberg_half)^(-elas)		


	minmat<-matrix(ncol=8,nrow=ncol(d2))
	
	for(i in 1:ncol(d2)){
    minmat[i,1]<-which.min(d2[,i])
    minmat[i,2]<-m2[minmat[i,1],i]
    minmat[i,3]<-s2[minmat[i,1],i]
    minmat[i,4]<-iceberg0[minmat[i,1],i]    
    minmat[i,5]<-which.max(mshares1[,i])   
    minmat[i,6]<-iceberg0[minmat[i,5],i]
    minmat[i,7]<-iceberg0[which.min(iceberg0[,i]),i]
    minmat[i,8]<-iceberg_half[which.min(iceberg0[,i]),i]
	}

	
	priceshock<-(1+minmat[,7])/(1+minmat[,8])
	phi0shock<-priceshock^(eps*(-0.1))
	

	summary(minmat[,2])
	summary(minmat[,3])
	summary(minmat[,4])
	summary(minmat[,6])
	sum(minmat[,4]==minmat[,6])/nrow(minmat)
	
	mean(colMins(iceberg0))
	median(colMins(iceberg0))
	mean(colMins(miceberg0))
	median(colMins(miceberg0))
	mean(colMins(riceberg0))
	median(colMins(riceberg0))
	
	m1_half<-dmat1*(tau_half)
	mshares1_half<-t(t(m1_half)/colSums(m1_half))
	v4$Phi_f_half<-as.numeric(colSums(m1_half))	
	
	m1_ruralhalf<-dmat1*(tau_ruralhalf)
	mshares1_ruralhalf<-t(t(m1_ruralhalf)/colSums(m1_ruralhalf))
	v4$Phi_f_ruralhalf<-as.numeric(colSums(m1_ruralhalf))
	
	m1_mainhalf<-dmat1*(tau_mainhalf)
	mshares1_mainhalf<-t(t(m1_mainhalf)/colSums(m1_mainhalf))
	v4$Phi_f_mainhalf<-as.numeric(colSums(m1_mainhalf))	

	
	#v4<-merge(v3,h,all.x=T)
	
	sv4<-v4[,c("region","district","ward","village_name","Phi_f_mainhalf","Phi_f_ruralhalf","Phi_f_half","Phi_f","weighted_dist_std","sampled")]
	
	sx<-read.dta13("Data/MaLT Calibration Replication.dta") 
	names(sx)[1]<-"region"
	sx$id<-seq(1,nrow(sx),1)
	
	sx<-merge(sx,sv4,all.x=T)

### Calculate Mark-ups

    elas0<-c(elas_IV)
    mult0<-c(2)
    
    count<-0	
    
    size<-18
    comps<-matrix(nrow=size*length(elas0)*length(mult0),ncol=11)	  
    
    for(elas in elas0){
	  for(mult in mult0){

		eps<-elas*mult
				
		dmatpre1<-as.matrix(init_d)
		dmatpre1<-dmatpre1/sum(dmatpre1)

		#Add deltas to columns
		dmat1<-matrix(dmatpre1,nrow=nrow(tau),ncol=ncol(tau))
		
		m1<-dmat1*(tau)
		mshares1<-t(t(m1)/colSums(m1))
		
		### Calculate composite shares for markup equation
		
		smatpre<-mshares1*Ebig
		smat<-smatpre/rowSums(smatpre)
		
		zi<-colSums(omega_mat*adopt_mat,na.rm=TRUE)
		
		Zmat<-matrix(NA,nrow=nrow(mshares1),ncol=ncol(mshares1))
		for(j in 1:nrow(Zmat)){Zmat[j,]<-zi}
		
		y2$mshare<-elas-elas*(eps-1)/eps*rowSums(smat*mshares1*(Zmat))
	
		y2$est_markup<-1/(y2$mshare)
		
		y2$mc<-y2$av_urea_rprice_USD_w1/(y2$est_markup+1)
		
		y2$Ts<-y2$deltas/y2$av_urea_rprice_USD_w1^(-elas)
		
		ones<-matrix(1,ncol=ncol(tau),nrow=1)		
	
		iceberg0<-tau^(-1/elas)-1
		iceberg_half<-iceberg0/2
		tau_half<-(1+iceberg_half)^(-elas)
		
		miceberg0<-mtau^(-1/elas)-1
		miceberg_half<-miceberg0/2
		riceberg0<-rtau^(-1/elas)-1
		riceberg_half<-riceberg0/2	
		
		tau_mainhalf<-(1+miceberg_half)^(-elas)*(1+riceberg0)^(-elas)	
		tau_ruralhalf<-(1+miceberg0)^(-elas)*(1+riceberg_half)^(-elas)		
		pricingfunc<-function(r){
			
			#r<-init_r
			#r<-y2$av_urea_rprice_USD_w1
			
			r<-ifelse(r==0,NA,r)
			
      #Calculate new delta's (quality adjusted prices) at new prices
			
			dmatpre<-y2$Ts*r^(-elas)
			dmatpre<-ifelse(is.na(dmatpre),0,dmatpre)
			
			### Calculate new lambdas and the shock

			dmat<-matrix(dmatpre,nrow=nrow(newtau),ncol=ncol(newtau))
			m<-dmat*(newtau)
			
			mshares<-t(t(m)/colSums(m))
			
			Shock<-(colSums(m)/Phi_f)/phi0shock2    ##### <<<<-------------- This is where we add-in the price shock
			
			### Calculate new adoption
			
			adopt_mat_new<-matrix(NA,nrow=nrow(adopt_mat),ncol=ncol(adopt_mat))
			for(k in 1:nrow(adopt_mat)){adopt_mat_new[k,]<-Shock/(Shock+(1-adopt_mat[k,])/adopt_mat[k,])}
			
			### Calculate new expenditures
			
			f2_new<-matrix(NA,nrow=nrow(f2),ncol=ncol(f2))
			for(k in 1:nrow(f2_new)){f2_new[k,]<-f2[k,]*Shock^(1/eps)*(adopt_mat_new[k,]/adopt_mat[k,])^((eps-1)/eps)}
			  
			omega_mat_new<-matrix(NA,nrow=nrow(omega_mat),ncol=ncol(omega_mat))
			for(i in 1:nrow(omega_mat_new)){
			  omega_mat_new[i,]<-f2_new[i,]/colSums(f2_new,na.rm=TRUE)
			} 

			### Calculate composite shares for markup equation

			zi_new<-colSums(omega_mat_new*adopt_mat_new,na.rm=TRUE)
			
			Zmat_new<-matrix(NA,nrow=nrow(mshares1),ncol=ncol(mshares1))
			for(j in 1:nrow(Zmat_new)){Zmat_new[j,]<-zi_new}			

			### Calculate new aggregate expenditures
			
			Ebig_new<-matrix(NA,nrow=nrow(Ebig),ncol=ncol(Ebig))
			
			for(j in 1:nrow(Ebig_new)){
			  Ebig_new[j,]<-colMeans(f2_new,na.rm=TRUE)*v3$census_bothsex
			}
			
			smatpre_new<-mshares*Ebig_new
			smat_new<-smatpre_new/rowSums(smatpre_new)			

			mshare_new<-elas-elas*(eps-1)/eps*rowSums(smat_new*mshares*(Zmat_new))
  		newmarkup<-1/(mshare_new)				

			X1<-(r-newmc*(1+newmarkup))
			X1<-ifelse(is.na(X1),0,X1)
			return(X1)		
			
		}

		outcomes<-function(r,cost,tn){

			#r<-ans$x
			#tn<-tau_half
			
			r<-ifelse(r==0,NA,r)
			
			#Calculate new delta's (quality adjusted prices) at new prices
			
			dmatpre<-y2$Ts*r^(-elas)
			dmatpre<-ifelse(is.na(dmatpre),0,dmatpre)
			
			### Calculate new lambdas and the shock
			
			dmat<-matrix(dmatpre,nrow=nrow(newtau),ncol=ncol(newtau))
			m<-dmat*(newtau)
			
			mshares<-t(t(m)/colSums(m))
			
			Shock<-(colSums(m)/Phi_f)/phi0shock2   ##### <<<<-------------- This is where we add-in the price shock
			
			### Calculate new adoption
			
			adopt_mat_new<-matrix(NA,nrow=nrow(adopt_mat),ncol=ncol(adopt_mat))
			for(k in 1:nrow(adopt_mat)){adopt_mat_new[k,]<-Shock/(Shock+(1-adopt_mat[k,])/adopt_mat[k,])}
			
			### Calculate new expenditures
			
			f2_new<-matrix(NA,nrow=nrow(f2),ncol=ncol(f2))
			for(k in 1:nrow(f2_new)){f2_new[k,]<-f2[k,]*Shock^(1/eps)*(adopt_mat_new[k,]/adopt_mat[k,])^((eps-1)/eps)}
			
			omega_mat_new<-matrix(NA,nrow=nrow(omega_mat),ncol=ncol(omega_mat))
			for(i in 1:nrow(omega_mat_new)){
			  omega_mat_new[i,]<-f2_new[i,]/colSums(f2_new,na.rm=TRUE)
			} 
			
			### Calculate new revenues
			
			Fi<-as.matrix(colSums(f2_new,na.rm=TRUE))

			### Calculate composite shares for markup equation
			
			zi_new<-colSums(omega_mat_new*adopt_mat_new,na.rm=TRUE)
			
			Zmat_new<-matrix(NA,nrow=nrow(mshares1),ncol=ncol(mshares1))
			for(j in 1:nrow(Zmat_new)){Zmat_new[j,]<-zi_new}			
			
			### Calculate new aggregate expenditures
			
			Ebig_new<-matrix(NA,nrow=nrow(Ebig),ncol=ncol(Ebig))
			
			for(j in 1:nrow(Ebig_new)){
			  Ebig_new[j,]<-colMeans(f2_new,na.rm=TRUE)*v3$census_bothsex
			}
			
			smatpre_new<-mshares*Ebig_new
			Rj<-rowSums(smatpre_new)			
			smat_new<-smatpre_new/rowSums(smatpre_new)			
					
			mshare_new<-elas-elas*(eps-1)/eps*rowSums(smat_new*mshares*(Zmat_new))
			marks<-1/(mshare_new)			
			
			newprofits<-Rj*(marks/(marks+1))

			return(list(adopt=(adopt_mat_new),expenditures=(f2_new),markup=as.numeric(marks),sales=as.numeric(Rj),profits=as.numeric(newprofits),Phihat=as.numeric(Shock)))

		}

		phi0shock2<-ifelse(phi0shock,1,1)   #### For baseline shock, set equal to one
		
		init_r<-ifelse(is.na(y2$av_urea_rprice_USD_w1),0,y2$av_urea_rprice_USD_w1)
		newmc<-y2$mc
		newtau<-tau
		ans<-nleqslv(init_r,fn=pricingfunc,control=list(allowSingular=TRUE,trace=1,maxit=10000))	
		adopt_mat_base<-outcomes(ans$x,y2$mc,newtau)$adopt
		exp_mat_base<-outcomes(ans$x,y2$mc,newtau)$expenditures
		rev_base<-outcomes(ans$x,y2$mc,newtau)$sales	
		profit_base<-outcomes(ans$x,y2$mc,newtau)$profits		
	  markup_base<-outcomes(ans$x,y2$mc,newtau)$markup
	  prices_base<-ans$x
	  shock_base<-outcomes(ans$x,y2$mc,newtau)$Phihat

	  phi0shock2<-phi0shock   #### Assign 10% price shock
	  
		init_r<-ifelse(is.na(y2$av_urea_rprice_USD_w1),0,y2$av_urea_rprice_USD_w1)
		newmc<-y2$mc
		newtau<-tau_half
		ans<-nleqslv(init_r,fn=pricingfunc,control=list(allowSingular=TRUE,trace=1,maxit=10000))	
		adopt_mat_half2<-outcomes(ans$x,y2$mc,newtau)$adopt
		exp_mat_half2<-outcomes(ans$x,y2$mc,newtau)$expenditures
		rev_half<-outcomes(ans$x,y2$mc,newtau)$sales
		profit_half<-outcomes(ans$x,y2$mc,newtau)$profits
		profit_half_noentry<-profit_half
		markup_half<-outcomes(ans$x,y2$mc,newtau)$markup
		prices_half<-ans$x
		shock_half<-outcomes(ans$x,y2$mc,newtau)$Phihat
		
		
		N<-50
		
		res<-matrix(nrow=N*nrow(v4),ncol=21)
		

		for(i in 1:nrow(v4)){
		
		  res[(1+(i-1)*N):(i*N),1]<-id[,i]
		  res[(1+(i-1)*N):(i*N),2]<-adopt_mat[,i]
		  res[(1+(i-1)*N):(i*N),3]<-f2[,i]
		  #res[(1+(i-1)*N):(i*N),3]<-v4$weighted_dist_std[i]
		  #res[(1+(i-1)*N):(i*N),4]<-v4$sampled[i]	  
		  res[(1+(i-1)*N):(i*N),14]<-i
		  res[(1+(i-1)*N):(i*N),4]<-adopt_mat_half2[,i]
		  res[(1+(i-1)*N):(i*N),5]<-exp_mat_half2[,i]	  
		  res[(1+(i-1)*N):(i*N),21]<-v4$weighted_dist_std[i]
		}
		
		res<-as.data.frame(res)
		names(res)<-c("id","adopt_base","exp_base","adopt_half","exp_half","adopt_halfmain","exp_halfmain","adopt_halfrural","exp_halfrural","adopt_halftrans","exp_halftrans","adopt_fisp","exp_fisp","village","shock_base","shock_half","shock_mainhalf","shock_ruralhalf","shock_halftrans","shock_fisp","weighted_dist_std")
		#names(res)<-c("id","adopt_base","exp_base","dist_std","sampled","adopt_half","exp_half")		
		
		sx2<-subset(res)

		base<-summary(lm_robust(log(adopt_base)~weighted_dist_std,cluster=village,sx2,se_type="stata"))		
		half<-summary(lm_robust(log(adopt_half/adopt_base)~weighted_dist_std,cluster=village,sx2,se_type="stata"))
		
		baseexp<-summary(lm_robust(log(exp_base)~weighted_dist_std,cluster=village,sx2,se_type="stata"))		
		halfexp<-summary(lm_robust(log(exp_half/exp_base)~weighted_dist_std,cluster=village,sx2,se_type="stata"))
	
		count<-count+1
		
	  }}
				

#Row 2 of appendix table B6
half
halfexp
    

    
